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We introduce a formulation of Eulerian general relativistic hydrodynamics which is applicable for 
(perfect) fluid data prescribed on either spacelike or null hypersurfaces. Simple explicit expressions 
for the characteristic speeds and fields are derived in the general case. A complete implementation 
of the formalism is developed in the case of spherical symmetry. The algorithm is tested in a number 
of different situations, predisposing for a range of possible applications. We consider the Riemann 
problem for a polytropic gas, with initial data given on a retarded/advanced time slice of Minkowski 
spacetime. We compute perfect fluid accretion onto a Schwarzschild black hole spacetime using 
ingoing null Eddington-Finkelstein coordinates. Tests of fluid evolution on dynamic background 
include constant density and TOV stars sliced along the radial null cones. Finally, we consider the 
accretion of self-gravitating matter onto a central black hole and the ensuing increase in the mass 
D ' of the black hole horizon. 
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I. INTRODUCTION 



The theoretical understanding of the detailed dynamics of black hole interactions with matter is today a key research 
activity, both in the gravitational wave 0] and high-energy astrophysics H communities. Such research targets the 
interpretation of data obtained (or anticipated) in a number of distinct observational windows, from high resolution 
qs ' spectra of regions suspected of harboring black holes, e.g., with satellite experiments like the Rossi X-ray Timing 
Explorer (RXTE) jy, to broad band gravitational wave detection efforts, e.g., the Laser Interferometric Gravitational 
■Xh Observatory (LIGO) . Computations within the Newtonian paradigm have reached high levels of sophistication (see, 
e.g., H), offering important clues and support for further development of astrophysical models. Clearly, computations 
within the framework of general relativity would be highly desirable. Several efforts are underway for meeting the 
difficult challenge of solving the Einstein equations in full generality (e.g., p-ROf). 

Yet the task before the computational scientists is daunting, as the underlying theory is rich in conceptual and 
technical complications. This attribute of relativistic gravity, led in the pre-60's period to a certain confusion with 
respect to the physical content of the theory, e.g., as to whether gravitational waves exist at all. The conundrum was 
finally resolved with the introduction and masterly manipulation of the characteristic initial value problem (CIVP) 
by Bondi and Sachs |ll],|l2|]. fmportant features of the CIVP suggest it may be serve as a valuable tool in the 
effort of studying interesting black hole spacetimes (a review of the development of algorithms based on the CIVP is 
given in [p.3|). We highlight, in particular, the long term stability of non-linear numerical evolution schemes based 
on null coordinates, shown for regular spacetimes in axisymmetry |L4J, and black hole spacetimes in full 3D |l5| . 
Additionally, the schemes are of low computational cost per spacetime point, an issue which acquires significance 
when stability problems are controlled. Both features ultimately derive from the gauge properties of the null foliation, 
which captures directly the wave degrees of freedom as a propagation equation for a complex function, with all other 
relevant equations converted into ODE's to be integrated along the characteristic surface. The CIVP formulation of 
the Einstein equations hence offers, with the inclusion of appropriate matter dynamics (e.g., in the form of perfect 
fluid hydrodynamics), among others, the possibility of accurate studies of generic black-hole-matter interactions. 

The aim of this paper is to present and test a general formulation of the general relativistic hydrodynamic (GRH) 
equations for ideal fluids, appropriate for numerical work, which is suited, but not restricted, to integration on 
spacetimes foliated with null hypersurfaces. The form invariance of the approach with respect to the nature of the 
foliation implies that existing work on highly specialized techniques for fluid dynamics can be adopted with minimal 
modifications. In our program of studying black hole matter interactions we have already used such techniques in 
estimating gravitational radiation from bulk fluid accretion J16| in the black hole perturbation limit. The developments 
here constitute a first step in extending that program into the fully non-linear regime. 



The paper contains four main sections. Section || presents the essential formal elements of a new prescription for 
solving GRH, i.e., a choice of conserved and primitive variables and their relationship, along with a diagonalization 



of the Jacobian matrix and explicit expressions for its eigenvalues and eigenvectors. In Section [II we formulate 
the problem of coupled evolution of matter and geometry fields in spherical symmetry. The analog of the Tolman- 
Oppcnhcimer-Volkoff (TOV) equation is constructed and integrated along the null cone. Those solutions are used 
for consistency and convergence tests. Those are demonstrated, along with some relevant details of the numerical 



implementation, in Section IV. That Section also discusses the important issue of shock-capturing, explored for fluid 
data posed on a null surface (in flat spacetime). In Section |V| we turn to black hole spacetimes and spherical accretion. 
The test fluid limit is used as a test-bed, and we proceed to a preliminary study of self gravitating accretion. 

We use geometrized units (G = c = 1) and the metric conventions of fll7|| . Boldface letters (and capital indices) 
denote vectors (and their components) in the fluid state space. 

II. FORMALISM FOR GENERAL RELATIVISTIC HYDRODYNAMICS 

A. Prelude to the formal developments 

The traditional approach for relativistic hydrodynamics on spacelike hypersurfaces is based on Wilson's pioneering 
work Jig ]. In this formalism the equations were originally written as a set of advection equations. This approach 
sidesteps an important guideline for the formulation of non-linear hyperbolic systems of equations, namely the preser- 
vation of their conservation form. This feature is necessary to guarantee correct evolution in regions of sharp entropy 
generation. Nevertheless, in the absence of such features, non-conservative formulations are equivalent to conservative 
ones. The approach is then simpler to implement and has been widely used over the years in a number of astrophysical 
scenarios. On the other hand, a numerical scheme in conservation form allows for shock-capturing, i.e., it guaran- 
tees the correct Rankine-Hugoniot (jump) conditions across discontinuities. Writing the relativistic hydrodynamic 
equations as a system of conservation laws, identifying the suitable vector of unknowns and building up an approxi- 
mate Riemann solver has permitted the extension of modern high-resolution shock- capturing (HRSC in the following) 



schemes from classical fluid dynamics into the realm of relativity |19|. The main theoretical ingredients to construct 



such a scheme in full general relativity can be found in |20l . An up-to-date collection of different applications of HRSC 



schemes in relativistic hydrodynamics is presented in plfl . We will return to these schemes in section IV below. 

Non finite-difference methods have also been applied recently to compute relativistic flows, most notably, pseudo- 
spectral methods and smoothed particle hydrodynamics (SPH) techniques. Although pseudo-spectral methods have 
enhanced accuracy in smooth regions of the solution, the correct modeling of discontinuous solutions is still their main 
drawback |22||. Recently, however, progress has been achieved in this direction by using multi-domain decomposition 
techniques |23|]. On the other hand, SPH, as any other particle method, suffers from being dissipative when resolving 
steep gradients |24|] . In spite of this, it has recently proven to work in the ultra-relativistic regime |g5[ . Additionally, it 
has been shown that it is possible to generalize SPH methods to hyperbolic systems other than the Euler equations |2q ] 

Procedures for integrating various forms of the hydrodynamic equations on null hypersurfaces have been presented 
before p^] (see [p8i for a recent implementation). This approach is geared towards smooth isentropic flows. A 
Lagrangian method, applicable in spherical symmetry, has been presented by p9[ . Recent work in BQ] includes an 
Eulerian, non-conservative, formulation for general fluids in null hypersurfaces and spherical symmetry, including 
their matching to a spacelike section. Here we show that the GRH equations can be formulated in a conservative and 
form invariant way (i.e., irrespective of the spacelike or null nature of the foliation) for an arbitrary three dimensional 
spacetime and any perfect fluid with polytropic equation of state. 

B. The formulation 

1. Variables and evolution equations 

We consider the relativistic conservation equations (continuity equation and Bianchi identities) upon introducing 
an explicit coordinate chart (x°,x l ), i.e., 

gJ» = 0, (1) 

gT^ = -^—gTl x T» x , (2) 




where the scalar x° represents a foliation of the spacetime with hypersurfaces (coordinatized by (x 1 , x 2 , x 3 )) , 

matter current and stress energy tensor for a perfect fluid given by J M = pu^", T^ v — phu^u v + pg^ v , where p is the 

pressure, u M is the fluid four velocity and h = I + e + p/p is the relativistic specific enthalpy. 

We do not restrict the foliation to be spacelike (that is, the level surfaces of x° may be also null). We define the 
coordinate components of the four- velocity w M = (u°,u % ). The velocity components w 4 , together with the rest-frame 
density and internal energy, p and e, provide a unique description of the state of the fluid and will be called (following 
common usage) the primitive variables. They constitute a vector in a five dimensional space w = (p,u l ,e). The 
index A is taken to run from zero to four, coinciding for the values (1,2,3) with the coordinate index i. We define the 
initial value problem for equations (0J2) in terms of another vector in the same fluid state space, namely the conserved 
variables, U B , individually denoted (D,S\E) f||], 

D = U° = J° = pu° , (3) 

S* = IP = T° l = phu°u l + pg m , (4) 

E = \J A = T 00 = phu°u Q +pg 00 . (5) 

With those definitions the equations will take the standard conservation law form [ p2[ , 

d x0 {J—gV A ) + a,, (V=gF' A ) = s A , (6) 

where \J—g is the volume element associated with the four-metric and we defined the flux vectors F J and the source 
terms S (which depend only on the metric, its derivatives and the undifferentiated stress energy tensor), 



(7) 



(8) 

The state of the fluid is uniquely described using either vector of variables, i.e., XJ A or xv A , and one can be obtained 
from the other, as will be shown later, via the definitions (§)-@ and the use of the normalization condition for the 
velocity vector, 

g^u v = -1 . (9) 

Specification of XJ A on the initial hypersurface, together with an equation of state (EOS) p — p{p, e), followed by 
a recovery of the primitive variables leads to the computation of the fluxes and source terms. Hence, the first time 
derivative of the data is obtained, which then leads to the formal propagation of the solution forward in time. No 
continuity of the data is required, since in practice the evolution is achieved with the (possibly approximate) solution 
of local Riemann problems. 

2. Local characteristic structure of the equations 

Utilizing the machinery of modern hydrodynamical methods, i.e., HRSC schemes, to integrate the previous equations 
requires the investigation of their local characteristic structure. We proceed here with this analysis. For this purpose 
we temporarily ignore the inhomogeneous part of the system. 

Introducing the Jacobian matrices 

B^ - ?*!± (10) 

the system can be written in quasi-linear form as 

^oU^+B^U^O. (11) 



F^o 


= J> = 


-pa 3 , 


F-J* 


_ rpji 


= phu l u° + pg tJ , 


FJ' 4 


— jy'O 


= phvPu 3 + pg° 3 


S° 


= o, 




S' 


= -V" 


~gTU T>lX > 


s 4 


= -yf- 


~9^lxT^. 



The characteristic speeds (eigenvalues) and eigenvectors of the system, in the j-th direction, are then given by the 
solution of the algebraic problems (we now omit the vector indices of the fluid space) 



det(W - X>1) = , 

(b j - xny = , 



(12) 
(13) 



where X J denotes the characteristic speed in the direction j, r J the corresponding eigenvector and I denotes the unit 
matrix. 

The strong coupling of the elements of the B matrix via the normalization condition hinders the analysis of the 
eigenvalue problem in terms of the conserved variables U. The customary method to proceed (see |3^| ), is to analyze 
the problem in terms of the primitive variables. Indeed, the suitable choice of the primitive variables, such that the 
eigenvalue problem may be solved explicitly, is of great practical value. 

With a choice of primitive variables w, the quasi-linear system (pT|) is rewritten as 

(14) 



A°3 a ,ow + A 3 d xJ w = , 



where 



A 



dU 
dw ' 

dw 



Hence, upon analyzing the algebraic system 



det{A J - P A ) = , 
{A? - P A V = , 



(15) 
(16) 



(17) 
(18) 



elementary algebra establishes that A- 7 ' = A- 7 and the eigenvectors of matrix B J are given by r J = A°r?'. 

We now proceed to diagonalize our system, with the choice of primitive variables w B = (p, u l ,e), and specializing 
the derivation to the case of a perfect fluid EOS, p — (T — l)pe, with T being the (constant) adiabatic index of the 
fluid. The procedure can be generalized to an arbitrary EOS and details will be given elsewhere. 

The matrices A are in this case 



A° = 



u° pfii 

hu°u k + (T - l)eg ok ph(vP5 k + u k n % ) pTu°u k + 



h(u ) 2 + (r - l)eg° 



(T-l)pg° k 
2phu°m pT(u ) 2 + (r - l)pg 00 



A 1 = 



u l pb) 

hu k u l + (T - l)eg kl phS k u l + phu k 8) pTu k u l + (T - l)pg kl 
hu°u z + (T - l)eg 0% phu Q 8) + phtfpj pTu°u l + (T - l)pg m 



with pi = ^j- = — ■^ i . For a given coordinate direction, which we label '1' (e.g., the x 1 direction), the matrix 
A 1 - A 1 A reads 



A 1 - A 1 A = 



a phi 

hu k a + (T - l)ec fc ph5 k a + phu k b % pTu k a + (T - \)pc k 
hu°a + (r — l)ed phvPbi + phapi pTvPa + (T — l)pd 



where the following shorthand notation is used 



a 



A 1 ^, 



..'■■i 



l„0fc 



yg 



5] - \ x m 



..m 



i„oo 



yg 



(19) 

(20) 



The eigenvalues of that matrix are 



X 1 - 



(triple) , 



and 



aV 



i 



Mc 2 s +v 1 (l - c 2 s ) T c sy /c 2 s M 2 + v l {\ - c 2 s ){2M - Cv^+Nil - cf(l - £)) 



± 1 " c 2 s (l - C) 
where all the metric information is encoded in the expressions C, M. and Af, 



„oo 



„01 



C 



K) 



(A 2 ' 



M 



(«°) 



cm • 



AA 



,ii 



K) 



0\2 ' 



Additionally, v 1 = \ and c s is the local sound speed satisfying 

hc 2 =X + -oK, 
P 2 

with x = 75^ = (r — l)e and k = ~2 = (r — l)p. A complete set of right-eigenvectors is given by 



dp 



>"o,i = u°(l,u 1 ,u 2 ,u 3 ,u°) . 



(21) 



(22) 



(23) 



(24) 



(25) 



r ,2 = (w° + p H i 3 2,u°u 1 + phu x p 32l u°u 2 + ph(u 2 p 32 - u°b 3 ), u°u 3 + ph(u 3 fi 32 + u°b 2 ), (u ) 2 + 2phu°/j 32 ) , (26) 



•"0,3 = {u° + pp 23 , u Q u l + phu x p 23n u°u 2 + ph{u 2 p 23 + u°b 3 ),u°u 3 + ph{u 3 p, 23 - u°b 2 ), (u ) 2 + 2phu°/j, 23 ) ■ 



(27) 



with the definitions 



r± = 



.£ u o B+ ff K 

e h 



P„0 



Bu l h + pY[A{u°(^ - ag 01 ) + u^K] 
-£u°Bu 2 h + pT[A{u a c 2 - ag 02 ) + u 2 K] 
~^u°Bu 3 h + pT[A(u°c 3 - ag 03 ) + u 3 K] 
-£(u°) 2 B + pT[2u°K - aAg 00 } 



A = 



u — u Hi 



d — c l ^ 
5 = 1 + 1, 



A = 



TaA 



r-i' 

B = l + hA, 
Hij = mb j - Hjh. 



(28) 

(29) 
(30) 



Note that the ± dichotomy in the last two eigenvectors is implicit in the corresponding non-degenerate (±) eigenvalues 
through the variables a, c l and d. 

The spectral decomposition given above applies to a chosen direction j. Since j is arbitrary, to obtain similar- 
expressions for the remaining directions, it suffices to specialize them accordingly, e.g., obtain the eigenvalues from 
expressions ( |2l|) and (|22]) with substitution of the desired direction, and permutation of the corresponding eigenvectors. 



3. Relations between variables 



So far the developments have been completely general. We specialize here the discussion to null coordinate systems. 

For the EOS commonly accepted, the propagation speeds of fluid signals are always sub-luminal. In addition, the 
bulk flow is always assumed to be a timelike vector field. Hence, the Cauchy initial value problem for the fluid is well 
defined for data given on a null hypersurface (see Fig. M. 



While the numerical algorithm updates the vector of conserved quantities (D, S Z ,E), we make extensive use of the 
primitive variables (p,u l ,s). Those would appear repeatedly in the solution procedure: in the characteristic fields, 
in the solution of the Riemann problem and in the computation of the numerical fluxes (see below). Hence, it is 
necessary to specify a procedure for recovering them from the conserved quantities. In the spacelike case the relation 
between the two sets of variables is implicit. An example of an iterative algorithm to recover the primitive variables 
in this situation can be found in |34|| . In the null case, the procedure of connecting primitive and conserved variables 
turns out to be explicit for a polytropic EOS. This is a direct consequence of the condition g 00 = which characterizes 
null foliations |12j| and leads to algebraic simplifications in the normalization expression g^u^u^ = — f . 

Assuming, hence, a perfect fluid EOS, the internal energy, e, can be directly obtained in terms of the conserved 
quantities as the positive solution to a binomial equation, more precisely 

A 2 

e = =, (31) 

D 2 + D V / D 2 +T(2-T)A 2 

where 

A 2 = -D 2 - g m E 2 - 2g 0l S l E - 9lJ S l S^ . (32) 

Once e is known the rest of the primitive variables follow, e.g., h = 1 + Te, p = D 2 h/E, and 

i S*-pg°* 



D(l+Te) 



(33) 



Having an explicit relation between conserved and primitive quantities has an impact on the efficiency of the 
numerical code, as it eliminates an iterative process that is required, at least once per each spacetime point. It is 
however unavoidable for general, e.g., tabulated, EOS. 

III. NON- VACUUM SPHERICALLY SYMMETRIC SPACETIMES IN THE TAMBURINO-WINICOUR 

FORMALISM 

We consider here the general spherically symmetric spacetime, with perfect fluid matter, following the formalism 



of Tamburino-Winicour 42 1, with a minor extension to cover the case of a worldtube immersed in matter. 



The Einstein equations sufficient for obtaining the spacetime development are grouped as 

(jr vr = Kl vr: W / 

{-J rv = fc-Lfrj v / 

G vv \r = nT vv \ r , (36) 

where the v coordinate is defined by the level surfaces of a null scalar (i.e., a scalar v satisfying V^vW^v = 0). 
The r coordinate is chosen to make the spheres of rotational symmetry have area 47rr 2 . The x 2 , x 3 coordinates 
in this geometry are simply taken to be the angular coordinates (#, <p) propagated along the generators of the null 
hypersurface, i.e., they parameterize the different light rays on the null cone. The first two equations contain only 
radial derivatives and are to be integrated along the null surface. The last equation (|36| ) is a conservation condition 
to be imposed on the world-tube V (see Fig. |l|). Equation ( p5| ) may be substituted for by the equivalent expression 
g ab R a b = 8Trg ab (T a b — g a i,T/2), where the indices (a, b) run over the remaining coordinates x 2 ,x 3 . To proceed with 
the integration, an additional choice of gauge must be made on the worldtube. This condition fixes the only remaining 
freedom in the coordinate system, namely the rate of flow of coordinate time at the world-tube. 

We present here the explicit expressions we use, including the matter terms. Adopting the Bondi-Sachs form of the 
metric element, 

e 2/3 V 

ds 2 = dv 2 + 2e 2l3 dvdr + r 2 {d9 2 + sinfl 2 ^ 2 ) , (37) 

r 

the geometry is completely described by the two functions /3(v,r) and V(v,r) (we will also interchangeably use the 
variable W = V — r). 

The (3 and V hypersurface equations are given by 

/3 r = 2nrT rr , (38) 

V r = e 2fi + AnrVT rr + 8nr 2 T rv , (39) 



the latter being equivalent (modulo the four-velocity normalization condition) to 

V r = e 2 f\l-A^r 2 {g AB T AB -T)). (40) 

The comma in the above equations indicates, as usual, partial differentiation. 

Boundary conditions (f3(v)r, V(v)r) for the radial integrations are provided by the equation (|36| ) which explicitly 
reads, 



V0,t , a Vv V r 1 



r 



2 



(Ve 










= 4:nrT vr , 

= 47rrT OT + 4Ttr 2 T vv /V , 

= 2wrT vr - 2Tre 2(3 T vv /g o 



, ,, £ + — e Ali - 4tt— T vv = , (41) 

r ^' 2F 2r 2r F 

with the adoption of a suitable gauge condition. For each choice one obtains a pair of ODE's along the worldtube. 
For example, with the choices 

V v = ~8nr 2 T vv , (42) 

(Ve 2 ? 
one obtains 

(43) 



respectively, where Vo and goo denote the integration constants at v = 0. All the above conditions are equivalent in 
vacuum, and since our present computations place the worldtube in very low density regions, we do not analyze the 
issue further. 

The hydrodynamic equations reduce, within our symmetry assumptions, to: 

D, v + F;° = -(lnV), v D - (lnV), r F r0 , (44) 

&v + Kr = -( ln V X^ r - (In V), r F rl - r^T"" , (45) 

E,v + Kr = -Q*V), V E - (\nV), r F ri - r^T"" , (46) 

where V = \J—g = r 2 sin^e 2 ^ is the four dimensional volume element. The precise form of the flux terms is obtained 
with direct use of the general formulae (0) and the Christoffel symbols V° nv are derived explicitly for metric (|37|). 

In summary, the initial value problem consists of equations (B8 39 M2-E6[) together with initial and boundary data 
for the fluid variables (p,e,u r ) on the initial slice So (at time vq) and the metric values of /3(uo)r ^md V^«o)r at the 
woldtube. Those equations and initial data are sufficient for obtaining the spacetime in a domain to the future of the 
initial hypersurface, which is radially bounded by the worldtube. 

A. Stationary configurations: Tolman-Oppenheimer-Volkoff solutions along the null cone 

We present here the equations describing a spherically symmetric equilibrium configuration in null coordinates. Such 
solutions constitute an excellent test-bed for consistency and accuracy checks of our algorithms. For the simplest 
derivation, it will be advantageous to use a slightly different form of the metric element. With the redefinition 
Y = Ve~ 2 P, it reads, 

e 4/3y 

ds 2 = dv 2 + 2e 2l3 dvdr + r 2 {dB 2 + sin 2 d<l> 2 ) . (47) 

r 

In analogy with the spacelike foliated case, the stationarity in all metric and fluid variables reduces the number of 
non-trivial equations to the following coupled pair, 

P,r=(^;-7^(l + 8irr 2 p)^ph 7 (48) 

Y r = l + 8Trr 2 {p-ph) 1 (49) 



where it is easily recognized that equation (fJ8|) is the radial momentum balance equation (|45|). This system of 
equations is solved as a system of ODE's. 

In the special case of a "stellar configuration", solutions are obtained by starting, e.g., with initial conditions 
p = p c ,Y = at the center of the star and integrating outwards until the pressure crosses the zero level [ [43[ . 
Following that recipe, the direct integration of Eq. (pq), with the boundary condition j3 = at the origin, completes 
the metric clement. Fields of a representative solution, which we use in the next section for code testing are presented 
in Fig. |. 

1. Constant density star 

In special cases the TOV equations can be integrated explicitly in terms of simple functions. A widely known 
example is the constant density star (po), derived here in ingoing null coordinates, 

/^m^, (50) 

V = ^-(3A-B), (51) 

B- A 

P = P0 3A—B> (52) 



where A = y/1 - 2M/R, B = y/l - 2Mr 2 /R 3 and M = 4irp R 3 , with R being the radius of the star. 

IV. ALGORITHMS AND TESTS 

The general formalism outlined in the previous two sections can form the basis for a variety of numerical approaches. 
Concerning the GRH equations, in the original Wilson's scheme [Eq], a combination of finite-difference upwind tech- 
niques with artificial viscosity terms were used to damp spurious oscillations, extending the classic treatment of shocks 
introduced by von Neumann (see, e.g., J3J|) into the relativistic regime. Artificial viscosity based methods, though, 
were later shown to fail at the threshold of the ultra-relativistic regime, once the Lorenz factor exceeds a value of 2. 
Explicit HRSC codes, following the so-called "Godunov approach" , appeared as a much more solid alternative. 

Since the early nineties, it has been gradually demonstrated (see, e.g., Q and references therein), that methods 
exploiting the hyperbolic character of the hydrodynamic equations are optimally suited for accurate integrations, 
even well inside the ultra-relativistic regime. As mentioned previously, these schemes are commonly known as high- 
resolution shock- capturing schemes. In a HRSC scheme, the knowledge of the characteristic fields (eigenvalues) of the 
equations, together with the corresponding eigenvectors, allows for accurate integrations, by means of either exact 
or approximate Riemann solvers, along the fluid characteristics. These solvers, which constitute the kernel of our 
numerical algorithm, compute, at every interface of the numerical grid, the solution of local Riemann problems (i.e., 
the simplest initial value problem with discontinuous initial data). Hence, HRSC schemes automatically guarantee 
that physical discontinuities appearing in the solution, e.g., shock waves, are treated consistently (the shock- capturing 
property). HRSC schemes are also known for giving stable and sharp discrete shock profiles. They have also a high 
order of accuracy, typically second order or more, in smooth parts of the solution. 

We proceed now to describe the highlights of our implementation. The grid structure is, in summary, as follows: 
An equidistant radial grid r^, with spacing Ar, denotes the location of cell centers on which the conserved variables 
and metric element components reside. The interface locations r\ = Ti — Ar/2 are used for the reconstruction of 
variables and the solution of the Riemann problems. Appropriate number of ghost- zones are added on the boundaries 
of the grid to allow imposition of boundary conditions. In particular, the worldtube conditions (J42||43|) are discretized 
in time along the first ghost-zone at rjy+i, where r^ denotes the last cell center inside the domain. The timestep is 
variable and consistently computed to satisfy the Courant condition for the hydrodynamic equations. The geometric 
equations, being ODE's in spherical symmetry, do not impose any timestep restriction. 

Our numerical implementation of the coupled system is broadly based on the notion of operator splitting, i.e., 
the integration of the initial value problem in steps, by successive applications of split components of the overall 
evolution operator. We describe here, schematically, the procedure. The set of equations, comprised of the three 
hydrodynamical equations in spherical symmetry, and the hypersurface equations governing the geometry, are written 
as 



a ( ,U + <9 r F(w,G) =0, (53) 

d v V =S[/(w,G), (54) 

A(U,w,G) = 0, (55) 

d r G =S G (w,G). (56) 

The ordering of the equations reflects the sequence in which they are solved in the code. We denote the metric variables 
collectively as G, while A stands for the set of algebraic relations connecting conserved and primitive variables (with 
the mediation of the metric). The velocity normalization condition is included with this set. 

The first step in the solution process involves an advection of the conserved variables U, from the initial data 
hypersurface So to the hypersurface Sau- The computation of the fluxes F uses the metric as it is given in So, along 
with the primitive variables that are assumed to have been computed on So hi the previous iteration. The advection 
may be performed by any modern numerical scheme for the propagation of non- linear waves, i.e., by any HRSC 
scheme. These schemes, being written in conservation form, are particularly well suited for this purpose. Hence, all 
conserved quantities in the differential equations are also conserved in their finite-differenced versions. More precisely, 
the time update of 

d v U + d r F = , (57) 

is done according to the following algorithm: 



Av 



ur +i = u;r- —( Ft+1/2 - Fl _ 1/2 ). (58) 

The index n represents the time level, while the time discretization interval is indicated by Av. The "hat" in the 
fluxes is used to denote the so-called numerical fluxes which, in a HRSC scheme, are computed according to some 
generic flux-formula, of the following functional form: 

F<±i = \ ( F ( U f±i) + F ( U f±i) " E I ^ I AS «M ' ( 59 ) 

Notice that the numerical flux is computed at cell interfaces (i ± 1/2). Indices L and R indicate the left and right 
sides of a given interface. The sum extends to p, the total number of equations. Finally, quantities A, Alo and r 
denote the eigenvalues, the jump of the characteristic variables and the eigenvectors, respectively, computed at the 
cell interfaces according to some suitable average of the state vector variables. 

In our code, the numerical integration of the hydrodynamic equations can be performed using two different ap- 
proximate Riemann solvers. These are the Roe solver J36J, widely employed in fluid dynamic simulations, with 
arithmetically averaged states (for the use of the Roe mean in a relativistic Roe solver see |37|) and the Marquina 
solver, recently proposed in |}8) (see also ]39J]). 

After the update of the transport terms the fluid variables are subsequently corrected for the effect of the source 
terms. Any stable ODE integrator is usable here. Our choice is a second order Runge-Kutta method. A point of 
interest here is that the source terms Sjj depend on both fluid and metric variables, and, in particular, on time 
derivatives of the latter. This implies that a second order capturing of the effect of those derivatives requires storing 
an additional time level of metric variables. 

The geometry equations comprise a system of radial ODE's, in which the right-hand-side depends on both the 
metric and the stress-energy tensor, a system to be solved on the hypersurface Sau- Initial conditions are provided on 
the world-tube. Hence, the integration cannot proceed without the recovery of the primitive variables on Sa-u, using 



the expressions given in section II B 3 . It can be seen that those expressions involve the metric (which is to be solved 
for). Furthermore, one must worry about the preservation of the velocity normalization condition, which strongly 
depends on the metric at a given point. This differential-algebraic entangling of metric and fluid variables, generic to 
any attempt to solve general fluid spacetimes using conservative formulations, is approached in the following way: 

The ODE integrator for equation (|5l) is chosen to belong to the implicit class. Our present choice is the second 
order Euler method. The metric variables in the new radial location r^+i are obtained iteratively (k denotes the 
iteration index). Inside the /c-iteration loop the intermediate metric variables G^ are used to obtain intermediate 
primitive variables w^ +1 , which then leads to updated values for the source term Sq i+1 . In addition, the zeroth 
component of the velocity vector is adjusted so that the primitive velocity field u M is normalized. This process can 
be described as 

A(U i+ i,wt lJ Gti) = J (60) 



A, 



G: + i-G, = -(S^ 1+ S GiI ), (61) 



2 



(the time level index has been suppressed here for clarity). In practice the number of required iterations is about four 
or five. The completion of this procedure furnishes simultaneously the metric and the primitive variables at the new 
radial location, at the new time level. 

A. The Riemann problem on a null Minkowski slice 

We proceed now to establish the feasibility of our proposal by performing numerical integrations of the GRH 
equations on null surfaces. Using the characteristic fields derived previously, we show here with two numerical 
examples the capabilities of existing HRSC schemes to integrate discontinuous solutions described in characteristic 
spacetime foliations. The initial data are given on a null slice (advanced or retarded time) of Minkowski space, and 
the evolution is compared to the suitably transformed exact solution. 

The shock tube (a particular case of the Riemann problem) is one of the standard tests to calibrate numerical 
schemes in classical fluid dynamics. The initial setup of this experiment consists of a zero velocity fluid having two 
different thcrmodynamical states on either side of an interface. When this interface is removed, the fluid evolves in 
such a way that four constant states occur. Each state is separated by one of three elementary waves: a shock wave, 
a contact discontinuity and a rarefaction wave. This time-dependent problem has an exact solution J4CJ], to which 
the numerical integration can be compared. In addition, it provides a severe test of the shock-capturing properties 
of any numerical scheme. In recent years, the shock tube problem began being used as a test of (special) relativistic 
hydrodynamical codes (see, e.g., p9| and references therein). The analytic solution of the Riemann problem, also 
available in relativistic hydrodynamics since |4l| , allows for a rapid and unambiguous comparison with the numerical 
evolutions. 

For our numerical demonstrations we consider two different initial setups. In case 1 the initial state of the fluid is 
specified by pl = 13.3, pi = 10 on the left side of the interface and pn = 0, pr = 1 on the right side. For numerical 
reasons, the pressure of the right state is set to a small finite value (pn = 0.66 • 10~ 6 ). Case 2 is as case 1 but with 
the left and right states reversed. We use a perfect fluid EOS with V = 5/3. 

1. The advanced time case 

In Fig. k| we plot the results for case 1. The left panels show the whole domain, with the x-coordinate ranging 
from to 1000. The right panels show a zoomed up view of the most interesting domain (x <G [750, 950]). We use a 
numerical grid of 2000 zones and, hence, a rather coarse spatial resolution (Ax = 0.5). 

^From top to bottom, Fig. [| displays the internal energy (e), the velocity (u x ) and the density (p). The thick 
dotted line represents the initial discontinuity. The solid line shows the exact solution after an advanced time v = 270 
and the dotted line indicates the numerical solution at this same time. In order to compute the exact solution in our 
new coordinates we have followed the procedure outlined in |41[] applying the appropriate time transformation, i.e., 
t = v — x. Clearly, the agreement between the exact (solid lines) and numeric solution is remarkable. The solution 
is characterized by an (outgoing) shock wave (moving to the right) and an (ingoing) rarefaction wave (moving to the 
left). One important point to notice is that features of the solution moving towards the left are more developed than 
those moving to the right. This is visible in the respective locations of the head of the rarefaction and the shock wave 
with respect to the initial discontinuity. The appearance of the rarefaction wave differs from the standard "spacelike" 
one (see, e.g., Fig. 2 of p9[]), due to the intervening time transformation. The internal energy and the density show 
an additional elementary wave, a contact discontinuity between the shock and rarefaction waves. Across the contact 
discontinuity pressure and velocity are constant, while the density exhibits a jump. 

In the right panels we note that the largest discrepancies are found at the tail of the rarefaction, in the form of 
an "undershooting" in density and internal energy and an "overshooting" in velocity. This is a typical feature of the 
numerical solution of the shock tube problem and it is not related to the null coordinates used here. We note that 
the constant state between the shock wave and the contact discontinuity is well resolved despite the coarse resolution 
used. For comparison purposes, the interested reader is referred to |39|| , where similar results were obtained (in a 
spacelike approach) employing a grid resolution 200 times finer. 

2. The retarded time case 

In Fig. [| the results for the "mirror" version of the shock tube problem 1 are plotted. As in Fig. g, the left panels 
show the whole domain, with the x-coordinate ranging from to 1000 whereas the right panels show a closed-up view 

10 



of the most interesting region (x S [180, 820]). We use the same grid resolution as in case 1. The initial data are now 
evolved up to a final time v = 120. Again, those features of the solution moving towards the left are more developed 
than those moving towards the right. The "tower" in the density is now considerably wider than in case 1 as well as 
the intermediate constant state in the internal energy between the shock and the contact discontinuity (almost not 
noticeable in case 1). This constant state is much more accurately resolved than in case 1 despite the use of the same 
resolution. The two corners of the rarefaction wave are also very well resolved. The major numerical deviations from 
the exact solution appear now at the post contact discontinuity region, with the density and the internal energy being 
slightly under and over estimated there. 

We should comment on the different ways in which the shock is captured in cases 1 and 2. Whereas in case 1 the 
shock is spread out in a number of grid zones (ss 10), in case 2 it is sharply resolved in 2-3 cells (out of the 2000 zones 
used). This is shown in Fig. pi. Those differences ultimately derive from the fact that a retarded (or advanced) time 
formulation of a Riemann problem (assuming, for the discussion, the initial discontinuity located at x — 0) breaks 
the commutation between the operation of reflection symmetry (between positive and negative values of x) and that 
of time evolution. The significance of this effect in applications involving shocks and its possible exploitation for 
ultrarelativistic outflows will be studied elsewhere. 

We also include in Fig. ra the evolution of the outgoing shock of case 1 at two different times, v = 300 and v = 600. 
For this run we are using 4000 zones with the x-coordinate ranging from to 2000. Hence, the spatial resolution is 
again Ax = 0.5. The initial discontinuity is placed at x = 1650. The solution was left to evolve for a time considerably 
longer than before in order to find out if the numerical representation of the shock had more time to accommodate 
itself into a sharper profile. From Fig. iwe find this not to be the case. The initial spreading present in the numerical 
solution remains constant in time. It is noteworthy to mention the independence of the number of points in which 
the shock is resolved from the total number of zones used. 

B. Convergence testing of the coupled code 

The TOV solutions presented before constitute a good test-bed for confirming both the consistency and the accuracy 
of the algorithm. The availability of semi-analytic solutions to the stationary problem means that there is a host of 
error functions that can be constructed by subtracting the evolved solution at a fixed total time tp from the initial 
profile. As an example, using the exact density profile pe, the quantity 

\\ P -Pe\\2 = Y.(p*-Pe) 2 (62) 

i 

measures the deviations produced by the numerical evolution. 

As our present implementations are geared towards black hole spacetimes (with topology of S 2 x R rather than 
regular R 3 spacetimes), in order to benefit from this test-bed without undue boundary complications, we consider 
only a portion of the spacetime, excluding the domain around r = 0. Initial data for the fluid and the metric are 
obtained with the solution of the TOV equations along the past null-cone of the center of symmetry. The integration 
then proceeds in the regular fashion described above, but restricted to a domain inside the star, e.g., in the sample 
shown in Fig. |2| this extends between r — 2 and r = 6. In the plots given in Fig. [7| we consider the convergence of 
such functions and their L^ norms to zero, as the grid spacing is reduced. The norm is found to converge to third 
order. This can be seen in the insert of Fig. 0, where the final values of the norm at v — 40 (about eight light crossing 
times) are plotted against grid size. This implies that the local error is second order. 

The constant density solution is obtained with the assumption of an ad-hoc (and largely unphysical) equation of 
state and hence it is not possible to evolve such data with the formalism developed here, but it is straightforward to 
test the integration of the hypersurface equations against the exact solution. We have confirmed that the integrated 
metric agrees along the hypersurface with the exact value to second order in the radial grid spacing. 

V. SPHERICAL ACCRETION ONTO A BLACK HOLE 

We proceed now in applying the formalism, in spherical symmetry, to the problem of interaction of matter with a 
black hole. In recent work [|44| the interaction of a scalar field with a black hole was investigated, partly as a probe 
into the concept of "singularity excision" [I45W47J. Following that same concept in spirit, in a previous investigation, 
we studied fluid interactions with the black hole geometry, in the test-fluid limit. We dubbed stationary coordinates 
that are regular at the horizon of an exact stationary black hole solution as horizon adapted coordinate systems pSy . In 
those coordinates the flow solution is smooth and regular at the black hole horizon. The steepness of the hydrodynamic 
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quantities dominates the solution only near the real singularity. This approach is now being applied in successively 
more general (and astrophysically interesting) fluid configurations [|49[|5C( | . We extend here this line of work by first 
performing test-fluid computations in the spirit of Jig] for the null Eddington-Finkelstein coordinate system. This is 
next further generalized to account for self-gravitating accretion flows. 

A. The test fluid limit 

We study spherical accretion of a perfect fluid onto a (static) black hole. The fluid is taken to have a sufficiently low 
density so that during the accretion process the mass of the black hole remains unchanged. Stationary solutions to 
this idealized problem can be computed exactly, up to algebraic equations. This was first derived for Newtonian flows 
by Bondi plj. The extension to general relativity was due to Michel p2|. The solution can easily be re-derived for 
coordinates other than the original Schwarzschild system employed by Michel (the details can be found in |48|| ), We 
use next this exact solution to quantify the accuracy of our numerical integrations. The significance of the test lies in 
its capturing of large curvature gradients near the black hole (i.e., it is a strong field computation), which translates 
into the existence of large source terms in the hydrodynamic equations. 

We use ingoing Eddington-Finkelstein (v,r,8,(f>) coordinates for which the line element reads 

ds 2 = - ( 1 - — J dv 2 + 2dvdr + r 2 (d8 2 + sin 2 Odcj) 2 ) , (63) 

where M is the mass of the black hole. The expressions for the characteristic fields of the hydrodynamic equations 
are then specialized using the above metric components. We use a polytropic EOS with Y (the adiabatic exponent 
of the gas) equal to 5/3. The numerical domain extends from any given non-zero radius inside the horizon, r m j„, to 
some outer radius r max (outside and far from the black hole horizon). In the particular simulation reported here we 
choose r m i n = 1.5M (inside the black hole horizon, located at r = 2M) and r max = 30M . We use a uniform (equally 
spaced) numerical grid of 200 zones. 

The test proceeds as follows: We set the semi-analytic solution as initial data throughout the domain, and then 
evolve these data, maintaining the exact solution as a constant inflow boundary condition at the outer boundary. 
Representative results are summarized in Figs. P~^0[ 

In Figure g we display only the solution up to a radius of 10M, in order to focus on the most interesting, strong 
field region, around the black hole horizon. The figure displays the primitive variables, (p, w r ,e), as functions of 
the ingoing Eddington-Finkelstein radial coordinate. These figures show the capturing of the steady-state spherical 
accretion solution. The solid lines represent the exact solution, while the filled circles indicate the numerical one. 
The latter has been evolved up to a time v = 500M. The agreement between the exact and numerical solutions is 
very good for all fields. This is more clearly seen in Fig. |9| where we plot the relative errors of the primitive variables. 
The largest errors appear at the innermost zones. The maximum error never exceeds 2%, for the internal energy, or 
1%, for the density, even with the relatively coarse grid of 200 points. The quality of the results for the velocity is 
excellent, with the maximum relative error being about 0.1%. 



In Fig. 10 the convergence properties are captured more quantitatively. The squared difference of the density from 
the exact result, integrated over the entire domain (the Li norm), is plotted as a function of advanced time, for a 
total time of 200. After an initial rise, the error settles into a final state, the level of which is converging to third 
order with the radial grid size. 

We note in passing that no secular instabilities of any kind arise during the computation, for total number of 
iterations of the order of 10 5 . Besides the intrinsic value of the test as an exact strong field solution, an important 
practical aspect merits mentioning here: very low density spherical inflow solutions constitute a good background flow, 
useful in computations where the primary dynamics of interest involves high density concentrations around the black 
hole. 

B. Accretion of a self-gravitating perfect fluid 

In our last numerical demonstration we investigate the accretion of a self-gravitating perfect fluid onto a non- 
rotating black hole. The setup of the simulation is as follows: Boundary values corresponding to a black hole of 
given mass Mp are specified at the world-tube, located at tb- The interior to that radius is filled first with low 
density data which do not interact with the black hole geometry in the timescales of interest. The dynamically 
interesting fluid component, typically a high density distribution (compared to the background) of compact support 
with sufficiently strong self-gravity, is added next. The data to be specified are (p, e, u r ) and are, in short, completely 
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arbitrary, with one exception: For the setup to correspond to a spacetime with a trapped region, i.e., with a horizon, 
the added fluid mass should not exceed Mp. The velocity profile is borrowed from the Bondi accretion solution and, 
hence, corresponds to an inwards monotonically increasing velocity. The initial data for the density are specified with 
explicit profiles, e.g., flat or Gaussian radial distributions. In our simulations we consider a Gaussian spherical shell 
surrounding the central black hole, with density parameterized according to 

P = Pb + Pme-^-^' (64) 

where pb is the background density. The rest of parameters take the values: p m = 10~ 4 , a = 0.1 and r c = 6M. 
The grid extends from r m i n — 1.1 Af to r max = 20 M. Finally, the internal energy is obtained assuming an initially 
isentropic distribution of pressure p = K p T (which would be valid at later times only for equilibrium configurations). 

The results of the simulation are plotted in Figs. O, O and H3. In Fig. O we display the evolution of the primitive 
variables (p, u r , e), from top to bottom, as a function of the ingoing Eddington-Finkelstein radial coordinate, r. The 
configuration is radially advected (accreted) towards the hole in the first 10M — 15M. Once the bulk of the accretion 
process ends, we are left with a quasi-stationary background solution (basically equivalent to the Bondi solution). In 
Fig. n3 we plot the evolution of the logarithm of the Riemann scalar curvature. We note that the initial shell has 
associated with it a non-negligible curvature. At late times the solution is again dominated by the curvature of the 
central black hole. 

The location of the apparent horizon of the black hole can be easily computed during the evolution. For the 
simplified case of spherical symmetry, this location is just given by the zero of the goo metric component. Our results 
show that the accretion process initiates a rapid increase of the mass of the apparent horizon. This is depicted in 
Fig. [13| The horizon almost doubles its size during the first 10M— 15M (this is enlarged in the insert of Fig. |l3|). Once 
the main accretion process has finished, the mass of the horizon slowly increases, in a quasi-steady manner, whose 
rate depends on the mass accretion rate imposed at the world-tube, T, of the integration domain. The numerical 
solution can be evolved as far as desired into the future. 

VI. SUMMARY AND CONCLUSIONS 

The equations of general relativistic hydrodynamics have been written directly in terms of components of the fluid 
fields in an arbitrary coordinate patch. The system of equations has been diagonalized in the general case, aiming at 
the subsequent use of advanced numerical methodology. The conservation form of the equations is preserved to the 
degree possible, whereas the formulation does not require a spacelike foliation for its implementation. 

We may ask whether there is any unexpected element in the developments here. The formulation of the hydro- 
dynamic equations as a Cauchy problem, (well) posed on a null spacetime surface, is intuitively expected and has 
been formulated before as mentioned in the text. The writing of the equations as a system of conservation laws on 
an arbitrary foliation is implicit in their abstract representation. It is further expected that for any non-degenerate 
choice of primitive fields, the system will be hyperbolic and diagonalizable. What does appear to be more of a happy 
algebraic coincidence rather than a general feature is the possibility of explicitly solving for the eigenfields. Indeed, 
choices of variables close to the ones presented here do not allow explicit diagonalization. We are aware of two cases 
in the literature where similar explicit resolutions as the ones reported here have been achieved. In the first case the 
authors made an explicit assumption of a spacelike foliation J2(J (see also II), while the second [pTJ appears moti- 
vated by a choice of variables specific to a non-relativistic Riemann solver (the Roe solver), and leads to expressions 
considerably more complicated than the ones presented here. 

We have developed a number of test-beds by recomputing known solutions in our coordinates. The performance of 
the algorithm was satisfactory in each instance, which establishes the overall feasibility of the approach. The practical 
value of this demonstration is that advanced schemes developed for the hydrodynamic equations by a large community 
of researchers (encompassing computational astrophysics) will be available with minimal modifications to this more 
specialized field. A further point of interest is that it is demonstrated here, implicitly, that the Newtonian pedigree 
of the field of hydrodynamics is essentially bypassed in the formulation of the relativistic equations as conservation 
laws. Our variables have no connection, in the null case, to the instantaneous rest frame (Eulerian) observers defined 
by the normal to a spacelike hypersurface. 

The selection of computations in the last section is geared towards highlighting the suitability of our approach to 
the study of black holes interacting with matter. In a representative computation, we have shown how a non-trivial 
increase of the mass of the black hole horizon can be achieved naturally in the present framework. This can be 
contrasted with the considerably harder task of achieving long term black hole evolutions in a spacelike approach, as 
it is seen, for example, in |53|. 



13 



Extending the present setup to three-dimensional spacetimes appears to be an important target for the near future. 
In this respect, the feasibility of extending the vacuum CIVP with the inclusion of matter sources has recently been 
reported p8| , 
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FIG. 1. Fluid equations of state are typically obeying causality, i.e., for equilibrium configurations the sound cone is contained 
within the light cone. In this case specifying fluid data on a null surface (here depicted as an advanced time null cone centered 
on the origin of coordinates) constitutes a Cauchy problem for the fluid. An example sound cone is depicted as a narrow 
forward cone, along with a worldline of a fluid element (in this case sub-sonic flow). This state of affairs would persist for an 
arbitrary curved spacetime and any null surface. 
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FIG. 2. The pressure and goi metric component for a representative TOV solution along the null cone (for T = 5/3, 
K — 4.349, central density p c = 8.1e — 4). A polytropic EOS is assumed, i.e., p = Kp T . 



16 



2.2 



1.4 



L. 
(1) 

c 
(1) 

i 

c 

(D 



■- 0.2 



-0.2 



12.( 
10.1 



1 6.( 
"° 4.1 



2.0 - 







- 






- 


- 


— initial data 
numeric 

— exact 








: / 


- 






- 


- 




- 




i 





1.2 

?« 
0.6 

| 0.4 

0.2 
+-> 

c 



1.2 -0.2 



- l.i 



0.6 i 


0.4 g. 

< 

0.2 



-0.2 6.0 



>4.0 

c 
(D 

"° 2.0 





1.2 



0.6 
0.4 
0.2 



< 






250.0 500.0 750.0 1000.1 

x-coordinate 



750.0 800.0 850.0 900.0 950.0 

x-coordinate 



FIG. 3. The outgoing shock tube problem (case 1). Exact versus numerical solution at an advanced time v = 270. The 
domain extends from x — to x = 1000 and a grid of 2000 zones is used (Aa; = 0.5). From top to bottom we plot the internal 
energy, velocity and density. The right panels show a zoomed view of the left panels focusing on the most interesting region. 
The thick dashed line shows the location of the initial discontinuity. The solid line is the exact solution and the dotted line is the 
numerically computed one. Those features of the solution moving to the right, i.e., the shock wave and contact discontinuity, 
are less developed than those moving to the left, i.e., the rarefaction wave. 
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FIG. 4. The ingoing shock tube problem (case 2). Exact versus numerical solution at an advanced time v = 120. The 
domain extends from x — to x — 1000 and a grid of 2000 zones is used (Ax = 0.5). From top to bottom we plot the internal 
energy, velocity and density. The right panels show a zoomed view of the left panels focusing on the most interesting region. 
The thick dashed line shows the location of the initial discontinuity. The solid line is the exact solution and the dotted line is 
the numerically computed one. As in Fig. H, those features of the solution moving to the right, i.e., the rarefaction wave, are 
less developed than those moving to the left, i.e., the shock wave and contact discontinuity. 
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FIG. 5. Capturing the shock wave. This figure shows the different ways the shock wave is resolved in the shock tube 
problems 1 (ingoing; left) and 2 (outgoing; right). In spite of the fact that the grid resolution is the same in both simulations, 
in case 2 the shock wave is spread out in a large number of cells whereas in case 1 it is sharply captured in two zones (out of a 
total number of 2000 zones). 
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FIG. 6. Location of the outgoing shock in the shock tube problem 1 at two different evolution times, v — 300 and v — 600. 
The numerical solution spreads the shock in a constant number of zones (sa 10) from the very start and throughout the 
evolution. 
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FIG. 7. The Li norm of the density function (i.e., the integrated squared difference of the numerical from the exact solution) 
is plotted against time, for four successive doublings of the grid size (curves A,B,C,D). The norm converges to third order. 
This can be seen in the insert, where the final values of the norm at v = 40 are plotted against grid size (diamonds). The best 
linear fit, represented by the line, has slope 3.05. This implies that the local error is second order. 
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FIG. 8. Perfect fluid spherical accretion in ingoing Eddington-Finkelstein coordinates (test fluid limit) . This figure compares 
the exact (solid lines) and numerical (filled circles) solutions for the hydrodynamical primitive variables, as a function of the 
radial coordinate. The numerical solution is evolved up to a time v = 500M. From top to bottom we plot the density, velocity 
and internal energy. The domain extends from 1.5M to 30M and a uniform grid of 200 zones is used. Only the first 10M are 
shown. 
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FIG. 9. Perfect fluid spherical accretion in ingoing Eddington-Finkelstein coordinates (test fluid limit): relative errors of 
the primitive variables. The maximum errors are less than 2% for the internal energy, 1% for the density and 0.1% for the 
velocity. These numbers correspond to a simulation using 200 radial zones in the interval 1.5M — 30M and for an evolution up 
tow = 500M into the future. 
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FIG. 10. Perfect fluid spherical accretion in ingoing Eddington-Finkelstein coordinates (test fluid limit): converging towards 
the stationary state. The four curves, labeled A,B,C,D, present the time evolution of the density L2 error norm (see Fig. M for 
definition), for successive doublings of the resolution. The insert shows the convergence of the norm at final time (v — 200Af ). 
The diamonds are measured values, the line is a linear best fit with slope 2.9. Of interest in this plot is the fast approach of 
the numerical solution to a steady state. 
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FIG. 11. Spherical accretion of a self-gravitating perfect fluid: evolution of the primitive variables. Ingoing Edding- 
ton-Finkelstein coordinates are used, in a grid of 500 zones spanning the radial interval between 1.1M and 20M. Three 
different times of the evolution are shown: v — (dotted thick line), v — 10M (solid thin line) and v = 15M (solid thick line). 
The spherical shell, centered at r — 6M, is radially advected towards the hole. A final quasi-stationary (Bondi) solution is 
achieved. 
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FIG. 12. Spherical accretion of a self-gravitating perfect fluid: evolution of the Riemann scalar curvature for the same 
simulation of Fig. O. Notice the non-vanishing self-gravity of the initial distribution, demonstrated by a curvature profile 
deviating significantly from the monotonic profile of a vacuum black hole. Once the shell is accreted, the solution is once again 
dominated by the curvature of the final vacuum black hole spacetime. 
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FIG. 13. Spherical accretion of a self-gravitating perfect fluid: evolution of the black hole apparent horizon mass. The mass 
of the apparent horizon shows a rapid increase in the first 15M (enlarged in the insert), in coincidence with the most dynamical 
accretion phase. The slow, quasi-steady growth at later times is the quiet response of the black hole to the low mass accretion 
rate imposed at the world tube. This rate can be made negligibly small (if desired). 
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